Example of Conditional Random Field Simulation Using the GMPE-SMTK and OpenQuake



In [ ]:

    
%load_ext autoreload
%autoreload 2
import warnings; warnings.filterwarnings("ignore")



In [ ]:

    
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm, Normalize
import smtk.hazard.conditional_simulation as csim
import smtk.sm_database_builder as sdb
from smtk.residuals.gmpe_residuals import Residuals
from smtk.residuals.residual_plotter import ResidualPlot, ResidualWithDistance
from smtk.parsers.sigma_database_parser import SigmaDatabaseMetadataReader, SigmaRecordParser, SigmaSpectraParser

Create Database from Event Records

record_folder = "data/LAquila_Good_Records/" database_folder = "data/LAquila_Database/" # Build the database builder = sdb.SMDatabaseBuilder(SigmaDatabaseMetadataReader, database_folder) builder.build_database("001", "LAquila Mainshock", record_folder) builder.parse_records(SigmaRecordParser, SigmaSpectraParser)

# Add on resolved horizonal components ims = ["PGA", "PGV", "Geometric"] sdb.add_horizontal_im(builder.database, ims)

Load in the database and retreive the GMPE residuals



In [ ]:

    
import cPickle
db1 = cPickle.load(open("data/LAquila_Database/metadatafile.pkl", "r"))



In [ ]:

    
# 1 GMPE, 2 IMS
gmpe_list = ["AkkarEtAlRjb2014"]
imts = ["PGA", "SA(0.2)", "SA(1.0)"]

resid1 = Residuals(gmpe_list, imts)
resid1.get_residuals(db1)



In [ ]:

    
ResidualWithDistance(resid1, "AkkarEtAlRjb2014", "PGA", plot_type="Linear", figure_size=(5,7))



In [ ]:

    
ResidualWithDistance(resid1, "AkkarEtAlRjb2014", "SA(1.0)", plot_type="Linear", figure_size=(5,7))

Load in the Rupture Model and View



In [ ]:

    
rupture_file = "data/laquila_rupture.xml"
rupture = csim.build_rupture_from_file(rupture_file)



In [ ]:

    
from openquake.hazardlib.geo.point import Point
from openquake.hazardlib.geo.polygon import Polygon
rupture_outline = []
for iloc in xrange(rupture.surface.mesh.shape[1]):
    rupture_outline.append(Point(rupture.surface.mesh.lons[0, iloc],
                                 rupture.surface.mesh.lats[0, iloc],
                                 rupture.surface.mesh.depths[0, iloc]))
for iloc in xrange(rupture.surface.mesh.shape[1]):
    rupture_outline.append(Point(rupture.surface.mesh.lons[-1, -(iloc + 1)],
                                 rupture.surface.mesh.lats[-1, -(iloc + 1)],
                                 rupture.surface.mesh.depths[-1, -(iloc + 1)]))
# Close the polygon
rupture_outline.append(Point(rupture.surface.mesh.lons[0, 0],
                             rupture.surface.mesh.lats[0, 0],
                             rupture.surface.mesh.depths[0, 0]))
rupture_outline = Polygon(rupture_outline)



In [ ]:

    
observed_sites = db1.get_site_collection()
pga_residuals = resid1.residuals["AkkarEtAlRjb2014"]["PGA"]["Intra event"]
#pga_residuals = (pga_residuals - (-3.0)) / 6.0
plt.figure(figsize=(10,8))
#ax = plt.subplot(111)
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(observed_sites.lons, observed_sites.lats, 
            s=40,
            c=pga_residuals, 
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.title("PGA Observed Intra-event Residual", fontsize=16)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)



In [ ]:

    
sa1_residuals = resid1.residuals["AkkarEtAlRjb2014"]["SA(1.0)"]["Intra event"]
#pga_residuals = (pga_residuals - (-3.0)) / 6.0
plt.figure(figsize=(10,8))
#ax = plt.subplot(111)
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(observed_sites.lons, observed_sites.lats, 
            s=40,
            c=sa1_residuals, 
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.title("Sa(1.0s) Observed Intra-event Residual", fontsize=16)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)



In [ ]:

    
# Generate a field of calculation sites
limits = [12.5, 15.0, 0.05, 40.5, 43.0, 0.05]
vs30 = 800.0
unknown_sites = csim.get_regular_site_collection(limits, vs30)

Generate a set of ground motion residuals conditioned upon the observations



In [ ]:

    
# Generate a set of residuals
output_resid = csim.conditional_simulation(observed_sites, pga_residuals, unknown_sites, "PGA", 1)



In [ ]:

    
plt.figure(figsize=(10,8))
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(unknown_sites.lons, unknown_sites.lats, 
            s=20,
            c=output_resid[:, 0].A.flatten(),
            marker="s",
            edgecolor="None",
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.scatter(observed_sites.lons, observed_sites.lats, 
            s=50,
            c=pga_residuals, 
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.title("PGA - Simulated Intra-event Residual", fontsize=16)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)



In [ ]:

    
sa1_residuals = resid1.residuals["AkkarEtAlRjb2014"]["SA(1.0)"]["Intra event"]
output_resid_1p0 = csim.conditional_simulation(observed_sites, sa1_residuals, unknown_sites, "SA(1.0)", 1)
plt.figure(figsize=(10,8))
#ax = plt.subplot(111)
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(unknown_sites.lons, unknown_sites.lats, 
            s=20,
            c=output_resid_1p0[:, 0].A.flatten(),
            marker="s",
            edgecolor="None",
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.scatter(observed_sites.lons, observed_sites.lats, 
            s=50,
            c=sa1_residuals, 
            norm=Normalize(vmin=-3.0, vmax=3.0))
plt.title("Sa (1.0s) - Simulated Intra-event Residual", fontsize=16)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)

Generate the Full Ground Motion Fields



In [ ]:

    
gmfs = csim.get_conditional_gmfs(db1,
                                 rupture, 
                                 sites=unknown_sites, 
                                 gsims=["AkkarEtAlRjb2014"],
                                 imts=["PGA", "SA(1.0)"],
                                 number_simulations=5,
                                 truncation_level=3.0)

Visualise the fields



In [ ]:

    
plt.figure(figsize=(10,8))

pga_field = gmfs["AkkarEtAlRjb2014"]["PGA"][:, 0]
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(unknown_sites.lons, unknown_sites.lats, 
            s=50,
            c=pga_field,
            marker="s",
            edgecolor="None",
            norm=LogNorm(vmin=0.001, vmax=1))
plt.xlim(12.5, 15.0)
plt.ylim(40.5, 43.0)
plt.title("PGA (g) - Conditional Random Field", fontsize=18)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)



In [ ]:

    
plt.figure(figsize=(10,8))
sa1_field = gmfs["AkkarEtAlRjb2014"]["SA(1.0)"][:, 0]
plt.plot(rupture_outline.lons, rupture_outline.lats, "r-")
plt.scatter(unknown_sites.lons, unknown_sites.lats, 
            s=50,
            c=sa1_field,
            marker="s",
            edgecolor="None",
            norm=LogNorm(vmin=0.001, vmax=1))
plt.xlim(12.5, 15.0)
plt.ylim(40.5, 43.0)
plt.title("Sa (1.0s) (g) - Conditional Random Field", fontsize=18)
plt.colorbar()
plt.xlabel("Longitude", fontsize=14)
plt.ylabel("Latitude", fontsize=14)